;;......................................................................................
;; NCL Graphics: 
;;......................................................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

begin

  load "load_setup.ncl"

  vname  = "CLDHGH" 
  vunits = "" 
  vscale = 1.0

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  casa  = (/"A: "+na, "B: "+nb/)
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)
  plotname = "figure_"+vname+"_map_"+shortname

  print(" ") 
  print(" ") 
  print("..................................................................") 
  print(" ") 
  print(" ") 
  print(" varname   : "+plotname) 
  print(" ctrl      : "+case(0)) 
  print(" test      : "+case(1)) 
  print(" plot name : "+plotname) 
  print(" ") 
  print(" ") 

  fngrid = "gridarea.nc"
  fgrid = addfile(fngrid,"r")
  gridarea = fgrid->cell_area

;;......................................................................................
;; Open netCDF file and read in dimensions
;;......................................................................................

   filename=case(0)
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))
   
;;......................................................................................
;; open wks 
;;......................................................................................
  wks   = gsn_open_wks (plotfmt,plotname)             ; open workstation


;;......................................................................................
;; colormap 
;;......................................................................................

  colormap = "BlWhRe"
  colormap = "WhiteBlueGreenYellowRed"
  colormap = "BlueYellowRed"
  colormap = "BlueYellowRedWhite"

  gsn_define_colormap(wks,colormap) 


;;......................................................................................
;; resources 
;;......................................................................................
  res1                      = True                ; plot mods desired
  res1@gsnFrame             = False               ; advance page
  res1@gsnDraw              = False                
  res1@gsnSpreadColors      = True                 ; use full colormap 
  res1@cnFillOn             = True                 ; turn on color
  res1@cnLinesOn            = False                ; turn on/off contour lines
  res1@cnFillMode           = "AreaFill"           ; default "AreaFill" "RasterFill"
  res1@lbLabelBarOn          = True                ; turn on/off the labelbar
  res1@lbOrientation         = "Vertical"         ; labelbar Orientation "Vertical""Horizontal" 
  res1@lbLabelsOn            = True                ; turn on/off the label of labelbar
  res1@lbBoxMinorExtentF     = 0.40                 ; area allotted to each box of the LabelBar    
  res1@lbLabelFontHeightF    = 0.02                                               
  res1@lbLabelAutoStride     = False                ; optimal labels         
  res1@cnLineLabelsOn        = False                ; turn on line labels
  res1@cnLinesOn             = False                ; turn on line 
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  ;res1@cnLevels             = (/0,10,20,30,40,50,60,70,80,90,100/) 
  res1@cnLevelFlags         = False
  res1@tiXAxisFontHeightF    = 0.020               ; XAxis Font
  res1@tiYAxisFontHeightF    = 0.018               ; YAxis Font
  res1@tiYAxisString         = " "
  res1@tiXAxisString         = " "
  res1@vpXF                  = 0.12     ; x location left
  res1@vpYF                  = 0.88     ; y location  up
  res1@vpWidthF              = 0.65     ; width
  res1@vpHeightF             = 0.35     ; height
  res1@gsnStringFontHeightF  = 0.02
  res1@tiMainString          = " "
  res1@tiMainFontHeightF     = 0.03
  res1@tiMainFontThicknessF  = 0.35
  ;res1@trYReverse            = True
  ;res1@gsnYAxisIrregular2Linear = True
  resP                 = True          
  resP@gsnPanelYWhiteSpacePercent = 2 ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 2  ; draw panel with white space added
  resP@txString            = vname
  resP@gsnPanelLabelBar    = False
  resP@lbOrientation       = "Horizontal" 
  ;resP@lbLabelFontHeightF  = 0.012               ; labels size
  ;resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 
  ;resP@lbLabelAngleF       = 15
  
;;......................................................................................
;; read in variables & create plot
;;......................................................................................

  plot = new(Ncase+1,graphic)
  
  do i=0,Ncase-1
     filename=case(i)
     fi=addfile(filename,"r")
    ;print(" read "+case(i))
     res1@gsnLeftString = casa(i)
     res1@gsnRightString= vunits                
     Var3dtema = fi->$vname$(time|:,lat|:, lon|:)
     Var3dtemp = Var3dtema(0,:,:) 
     Var3dtemp=Var3dtemp*vscale
     global_total = sum(Var3dtemp*gridarea)/sum(gridarea) 
     res1@gsnCenterString = sprintf("%5.3f", global_total) 
     plot(i)= gsn_csm_contour_map_ce(wks,Var3dtemp(lat|:,lon|:), res1)  
     if(i.eq.0) then
        vara = Var3dtemp
        aa = global_total
     end if
     if(i.eq.1) then
        varb = Var3dtemp
        bb = global_total
     end if
     delete(Var3dtemp)
  end do

  varc = vara
  varc = varb - vara

  delete(res1@cnLevels)

  res1@gsnLeftString = "B-A"
  res1@cnLevels = (/-0.1,-0.08,-0.06,-0.04,-0.02,0.02,0.04,0.06,0.08,0.1/) 
  global_total = sum(varc*gridarea)/sum(gridarea)
  res1@gsnCenterString = sprintf("%5.3f", global_total) 


  plot(2) = gsn_csm_contour_map_ce(wks,varc(lat|:,lon|:), res1)
  
  gsn_panel(wks,plot,(/3,1/),resP)
  
end


